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ABSTRACT 

The  purpose  of  this  study  is  to  demonstrate  the  feasibility  of  modeling  inelastic 
responses  of  ceramic  matrix  composites  (CMC),  where  deviations  from  elastic 
behavior  are  primarily  due  to  micro-crack  evolution  within  the  material.  The 
theoretical  basis  for  this  research  is  a  physically  based  homogenization  theory,  termed 
“discrete  damage  space  homogenization  method”  (DDSHM),  which  combines 
micromechanics  and  thermodynamics  to  determine  the  overall  response  functions  of 
multi-phase  materials  of  arbitrary  complexity.  The  procedure  is  named  DDSHM 
because  the  constitutive  and  evolution  equations  of  the  composite  are  obtained  for 
discrete  values  of  damage  (numbers  of  cracks  and  crack  lengths)  and  not  in  terms  of 
an  average  damage  parameter.  The  micro-crack  evolution  law  used  may  in  general 
include  time  and  temperature  effects  as  well  as  other  material  parameters  including 
changes  in  constituent  and  interphase  strengths  under  extended  environmental 
exposures. 

In  this  methodology,  a  thermodynamically  based  constitutive  relation  is 
formulated  using  two  scalar  functions:  a  thermodynamic  potential  which  specifies  the 
state  of  the  material  point;  and  a  dissipation  potential  which  governs  the  evolution. 
Given  the  constitutive  relations  of  the  individual  phases,  both  the  thermodynamic  and 
dissipation  potentials  are  completely  derived  by  first  solving  a  micromechanics 
(boundary  value)  problem,  and  then  performing  a  homogenization  procedure  which 
averages  the  solution  over  the  chosen  RVE.  The  constitutive  and  evolution  equations 
(CEE)  for  the  materials  to  be  studied  are  taken  to  be  functions  of  a  set  of  internal  state 
variables  and  a  set  of  damage  parameters  whose  number  must  be  selected  based  upon 
some  knowledge  of  the  behavior  of  the  micro-structure  under  load.  The  derived  CEE 
form  a  system  of  ordinary  differential  equations  that  describe  the  evolution  of  the 
micro-structure  of  the  RVE.  As  a  result,  the  FEM  based  micromechanics  analysis  can 
be  performed  independently  of  the  global  structural  analysis.  The  overall  approach  for 
deriving  the  constitutive  and  evolution  equations  has  already  been  integrated  with  the 
commercial  finite  element  code  ABAQUS  for  analysis  of  layered  polymer  matrix 
composites  and  is  extended  here  for  application  to  CMCs. 


C.  A.  Meyers,  A.  A.  Caiazzo,  J.  Meeker,  Materials  Sciences  Corp.,  135  Rock  Road,  Horsham, 
Pennsylvania  19044,  U.S.A. 
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INTRODUCTION 


Ceramic  matrix  composites  (CMCs)  represent  an  important  enabling  technology 
for  advanced  aerospace  engine  components.  Significant  progress  has  been  made  in 
understanding  the  behavior  of  these  complex  materials  in  the  linear  elastic  range,  and 
with  some  simple  forms  of  progressing  damage.  However,  a  general  physics-based 
methodology  to  assess  the  initial  and  lifetime  response  of  complex  CMC  material 
architectures  with  evolving  damage  would  provide  a  step  forward  in  terms  of 
acceptance  and  implementation  of  CMCs  in  aerospace  applications. 

For  the  purpose  of  discussion  we  assume  that  the  global  analysis  of  the 
composite  structure  is  to  be  conducted  using  the  Finite  Element  Method  (FEM). 
References  [1-4]  present  a  homogenization  theory  and  computational  strategy  that 
combines  micromechanics  and  thermodynamics  to  allow  determination  of  the 
overall  response  functions  of  multi-phase  materials  of  arbitrary  complexity.  In  this 
work,  a  thermodynamically  based  constitutive  relation  is  formulated  using  two 
scalar  functions:  a  thermodynamic  potential  that  specifies  the  state  of  the  material 
point;  and  a  dissipation  potential  which  governs  the  evolution.  Given  the 
constitutive  relations  of  the  individual  phases,  both  the  thermodynamic  and 
dissipation  potentials  can  be  completely  derived  by  first  solving  a  micromechanics 
problem,  and  then  performing  a  homogenization  procedure  which  averages  the 
solution  over  the  chosen  representative  volume  element  (RVE).  Even  in  the 
presence  of  growing  damages,  the  effective  constitutive  equations  need  be 
determined  only  once  and  independently  of  the  thermo-mechanical  load  path  that 
will  later  be  imposed  during  structural  analysis  calculations.  In  other  words,  this 
work  presents  a  computational  strategy  for  the  derivation  of  constitutive  equations 
for  composites  with  specific  growing  cracks  that  are  obtained  once  and  that  are 
valid  for  any  load  history  while  maintaining  consistency  with  continuum 
thermodynamics.  The  approach  has  been  termed  a  discrete  damage  space 
homogenization  method  (DDSHM)  because  the  effective  constitutive  and  evolution 
equations  (CEE)  of  the  composite  are  obtained  for  discrete  values  of  damage 
(numbers  of  cracks  and  crack  lengths)  not  in  terms  of  an  average  damage  parameter. 

The  goal  of  DDSHM  is  to  deliver  an  efficient  computational  procedure  for  the 
analysis  of  composite  materials  with  growing  cracks  that  retains  some  of  the 
benefits  of  homogenization  theory  and  micromechanics.  In  particular,  a  strategy  is 
sought  to  derive  the  fracture  mechanics  based  CEE  for  a  given  RVE  of  material. 
The  data  for  these  equations  are  computed  independently  from  the  global 
calculation,  i.e.,  the  detailed  RVE  problem  is  not  embedded  in  the  global  model. 
The  results  are  valid  for  any  load  history  and,  therefore,  the  constitutive  equations 
can  be  used  in  any  global  analysis.  This  allows  the  nonlinear  material  model  to  be 
integrated  into  a  commercial  finite  element  code  for  the  structural  analysis.  This 
has  been  accomplished  using  the  general  purpose  structural  analysis  code  ABAQUS, 
via  user-material  subroutines  (UMAT)  [5]. 
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THEORY 


In  this  section  of  the  paper,  the  theory  of  the  DDSHM  is  briefly  discussed  as 
pertaining  to  the  macroscopic  failure  of  composites.  DDSHM  is  a  computational 
strategy  used  to  solve  the  constitutive  and  evolution  equations  for  a  linearly  elastic 
composite  with  evolving  micro-cracks.  This  physically  based  method  stems  from 
the  homogenization  theory  presented  in  [1-4]  and  fracture  mechanics.  In  particular, 
the  effective  behavior  of  the  composite  material  is  governed  by  two  scalar  functions: 
a  thermodynamic  potential  and  a  dissipation  potential.  The  thermodynamic 
potential,  namely  the  Helmohotz  free  energy  H,  is  a  function  of  the  selected 
macroscopic  state  variables: 

H(t)  =  H(Eij(t),Zi(t),-,ZN(t))  (1) 


where  Ey(t)  is  the  small  strain  tensor  and  A|(t),...,  ^N(t)  are  N  internal  state 
variables  (IS Vs),  and  t  is  time.  The  actual  constitutive  equations  are  given  by 


v  dH  dH 

l,..  = -  and  (j,  = - 

9  8E.  1  dZ 


(2) 


where  2;j  is  the  macroscopic  stress  tensor  and  G;  is  the  generalized  thermodynamic 
force  promoting  the  growth  of  the  i,h  ISV.  The  dissipation  potential  Q  is  a  function 
of  the  generalized  thermodynamic  forces  G,: 


Q  =  Q(Gf),  i  =  , 

and  the  actual  evolution  equation  for  the  ith  ISV  is  given  by 


(3) 


Z  = 


da 

8G: 


(4) 


Consider  a  linearly  elastic  composite  with  growing  micro-cracks,  the  IS  Vs  can 
be  chosen  as  the  sizes  of  the  N  cracks  present  in  the  selected  representative  volume 
element  (RVE);  hence,  in  this  case,  G;  represents  the  energy  release  rate  for  the  i,h 
crack.  From  the  principles  of  linear  elastic  fracture  mechanics,  the  dissipation 
potential  can  be  expressed  as  the  following  form: 


n  =  \’l,{G,-G?)\i  =  l,-,N,  (5) 

where  q,  and  G,cr  are  referred  to  as  the  crack  growth  viscosity  coefficient  and  the 
critical  energy  release  rate  for  the  ith  crack,  respectively,  and  (•)  denotes  the  positive 
part  operator  defined  as 


|0,  if(f><  0 

\<f>,  if<f>>0’ 


(/)  e  91 . 


(6) 
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Moreover,  Eq.  (5)  yields  a  crack  evolution  law  for  the  /th  crack: 


i,=7,(G,-Gr>.  (7) 

In  order  to  use  the  above  theory,  H  must  be  evaluated  for  all  possible  values  of  the 
state  variables.  In  general,  it  is  impractical  if  not  impossible  to  accomplish  this  task  in 
an  exact  sense.  In  the  DDHSM,  H  is  evaluated  at  a  discrete  number  of  values  of  the 
ISVs,  and  interpolation  is  then  carried  out  to  obtain  the  general  function,  e.g.  Lagrange 
polynomials.  Specifically,  for  a  given  crack  state,  i.e.  the  values  of  X,  at  time  t,  and 
based  on  the  effective  constitutive  equation  for  the  composite  Eq.  (8), 

(t)  =  Cpqrs(\{t),X2(t),---,  XN  (. t))E„  (0 ,  (8) 

H  can  be  expressed  in  terms  of  the  effective  moduli  and  macroscopic  strain  as: 

Hit)  (9) 

Here,  Cyu  is  the  effective  stiffness  matrix  for  the  composite,  which  is  evaluated 
from  micromechanics  solution  for  a  representative  volume  element  (RVE). 

From  Eq.  (2)  and  (7),  a  system  of  ordinary  differential  equations  that  describe 
the  RVE  microstructure  evolution  in  terms  of  the  macroscopic  applied  strain  history 
can  be  formed;  the  evolution  of  the  independent  state  variabless  at  any  point  in  a 
loading  history  and  the  associated  overall  material  properties  can  then  be 
determined,  given  the  critical  energy  release  rates.  In  addition,  the  effective 
constitutive  equations  are  determined  only  once  and  independent  of  the  thermo¬ 
mechanical  load  path;  given  an  arbitrary  loading  condition,  the  associated  overall 
material  properties  can  also  be  determined.  It  is  also  noted  that  the  effective 
constitutive  equations  need  to  be  evaluated  at  various  stages  of  crack  development; 
hence,  all  possible  crack  paths  within  the  chosen  RVE  must  be  anticipated  a  priori. 


MODELING  APPROACH 

Three  basic  steps  are  required  for  any  constitutive  modeling  approach:  (i)  a  pre¬ 
processing  phase,  where  a  model  for  the  local  material  behavior  is  created;  (ii)  a 
solution  phase,  where  the  global  structural  behavior  is  determined  using 
information  on  the  evolving  material  state  from  the  constitutive  law;  and  (iii)  a 
post-processing  phase,  where  material  state  data  is  stored  for  review. 

The  DDSHM  can  be  viewed  as  a  hybrid  methodology  which  combines  classical 
homogenization  theory  with  the  principles  of  global-local  analyses  to  numerically 
derive  effective  constitutive  equations  that  are  load  history  independent  and  can  be 
obtained  separately  from  a  specific  global  analysis.  Indeed,  this  procedure  may  be 
viewed  as  a  computationally  efficient  alternative  to  traditional  global-local 
approaches  for  modeling  the  effective  behavior  of  composite  materials  with 
growing  cracks.  Items  1-3  (noted  by  boxes)  in  Figure  1  highlight  the  differences 
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between  the  DDSHM  and  the  global-local  approach.  Item  1  is  a  pre-processing  step 
which  is  required  by  the  DDSHM  but  not  by  the  global-local  method.  Note  that  the 
fact  that  this  extra  step  is  performed  in  the  preprocessing  phase,  i.e.,  decoupled 
from  the  global  solution  of  later  structural  analyses,  is  perhaps  the  most  unique  and 
desirable  aspect  of  the  DDSHM.  The  parameter  nSOL  depends  on  the  damage 
space  discretization,  i.e.,  combines  the  number  of  possible  crack  paths  as  well  as  the 
number  of  points  along  each  path  at  which  the  effective  elastic  moduli  are  explicitly 
calculated.  As  described  previously,  this  step  is  currently  accomplished  using  a 
finite  element  program  written  especially  for  performing  the  DDSHM  procedure. 

The  major  computational  difference  between  the  two  methods  occurs  during  the 
solution  phase  (Items  2  and  3).  Crack  evolution  in  DDSHM  is  determined  by 
numerically  solving  a  set  of  N  ordinary  differential  equations  (ODE’s)  indicated  by 
the  Griffith-type  crack  evolution  law.  This  is  currently  accomplished  using  a  4th 
order  Runge-Kutta  algorithm  with  adaptive  (time)  step  size  control.  In  the  global- 
local  method,  crack  evolution  is  determined  by  performing  a  complete 
micromechanics  analysis  for  each  material  point  in  the  global  analysis.  Hence,  if  an 
identical  RVE  is  selected  for  both  the  DDSHM  and  a  global-local  analysis,  the 
comparison  of  interest  is  the  computational  time  required  for  numerical  solution  of 
N  ODEs  versus  solving  a  mDOF  set  of  linear  equations,  where  mDOF  represents 
the  number  of  degrees  of  freedom  (DOF)  in  the  micro-mechanical  solution  of  the 
RVE  problem.  Solving  a  system  of  N  coupled  first  order  differential  equations, 
where  in  practice  N  will  in  general  be  on  the  order  10,  will  always  require 
significantly  less  time  and  computational  effort  than  that  needed  to  carry  out  the 
finite  element  analysis  of  size  mDOF,  even  if  the  micromechanics  finite  element 
analysis  can  be  conducted  without  user  intervention1.  In  other  words,  for  practical 
problems  the  computational  time  associated  with  performing  a  micromechanics 
analysis  which  involves  carrying  out  a  solution  for  a  FE  model  p  *  q  x  mDOF  times 
(i.e.,  for  each  global  load  increment)  is  expected  to  be  so  large  that  it  will 
considerably  outweigh  the  DDSHM  requirement  of  having  to  solve  nSOF  problems 
to  first  generate  the  set  of  N  ODEs  which  is  then  solved  p  x  q  times.  Furthermore, 
once  the  CEE  are  determined  by  the  DDSHM  they  can  be  applied  to  any  global 
analysis  problem  of  interest.  That  is,  the  DDSHM  computations  carried  out  in  the 
pre-processing  phase  are  not  unique  to  the  global  analysis  to  be  conducted.  Results 
of  a  global-local  analysis  cannot  be  re-used  for  any  other  problem  than  that 
explicitly  solved. 


1  As  reported  in  the  literature,  global-local  analyses  are  commonly  conducted  by  literally  extracting 
force  or  displacement  results  from  a  global  analysis  and  manually  applying  them  to  a  local  model  to 
assess  failure.  However,  it  is  possible  to  automate  the  process.  In  this  case  the  UMAT  routine  would  be 
a  FE  program. 
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Damage  Modeling  Approach 
DSSHM  Global-Local 


Figure  1 .  A  flow  chart  showing  computational  steps  for  a  nonlinear  structural 
analysis.  Differences  between  the  DDSHM  (left)  and  a  global-local  approach  (right) 
are  highlighted.  The  steps  shown  in  the  center  are  common  in  the  two  procedures. 

The  DDSHM  approach  automatically  yields  the  energy  available  to  grow  each 
possible  crack  for  each  material  point  in  the  analysis  based  upon  the  strain  state  acting 
at  that  point.  This  is  very  different  than  the  global-local  approach  where  results  of  a 
global  stress  analysis  are  reviewed  to  determine  the  areas  where  a  local  VCCT 
analysis  is  to  be  conducted.  If  the  material  properties  (Gcr  and  elastic  moduli)  are 
assumed  to  be  spatially  independent,  results  using  the  DDSHM  approach  immediately 
indicate  where  damage  onset  is  expected.  The  key  to  the  computational  efficiency  of 
the  DDSHM  is  the  manner  in  which  the  available  energy  release  rate  G;  is  practically 
computed.  The  DDSHM  derived  CEE  reduce  to  a  system  of  ordinary  differential 
equations  that  describe  the  RVE  microstructure  evolution  in  terms  of  the  macroscopic 
applied  average  strain  history.  Thus,  the  value  of  the  IS  Vs  at  any  point  in  a  loading 
program  can  be  determined  through  numerical  solution  of  the  set  of  ODEs  and  not  by 
repeatedly  solving  a  micromechanics  BVP  for  each  load  state.  The  DDSHM  is 
general  and  subject  only  to  those  limitations  associated  with  using  a  RVE  approach  for 
determination  of  effective  material  properties. 

Because  the  first  step  in  this  approach  is  to  build  and  run  a  RVE  with 
appropriate  material  properties  and  representations  of  initial  crack  surfaces,  one 
must  determine  whether  or  not  the  number  and  location  of  patches  and  available 
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crack  paths  is  sufficiently  dense  to  predict  damage  evolution  in  the  RYE  for  all 
possible  load  histories.  Microstructural  data  on  initial  configuration  and  damage 
progression  in  CMCs  were  examined  in  order  to  rationally  define  cracks  paths  and 
identify  any  CMC  characteristics  that  would  have  to  be  incorporated  in  the  code 
already  developed  for  PMCs.  In  particular,  micrographs  for  S200-1  which  had  been 
selected  as  the  focus  for  this  study  [6,  7],  indicated  that  the  primary  differences 
were  in  greater  porosity  and  larger  interface/interphase.  Images  of  microcracking 
and  fiber  pull-out  were  similar  (c.f.  Figure  2)  indicating  that  modeling  strategies 
used  previously  for  GRP  might  be  applied  with  appropriate  modification.  In 
particular,  an  existing  hexagonal  array  was  adapted  to  include  a  BN  interphase  and 
potential  sleeve-like  cracks  running  along  the  matrix/interphase  and  fiber/interphase 
boundaries. 

In  order  to  facilitate  the  creation  and  meshing  of  RVEs,  Materials  Sciences  Corp. 
has  developed  an  add-on  for  ABAQUS  CAE.  The  add-on  was  necessary  since  most 
FE  pre-processors  cannot  cleanly  mesh  around  internal  delaminations.  Using  this  add¬ 
on,  a  delamination  may  be  any  open,  non-intersecting  surface  which  the  user  defines 
by  a  series  of  faceted  open  curves.  The  add-on  lofts  straight  lines  between  the  vertices 
of  the  curves  to  form  the  edges  of  the  facets.  The  add-on  then  builds  a  geometric 
representation  of  the  RVE  and  assigns  mesh  seeds  to  the  crack  surfaces  so  that  surface 
nodes  will  only  be  placed  at  the  facet  vertices.  If  desired,  the  user  can  then  intervene 
and  manually  use  ABAQUS  CAE’s  meshing  tools  to  build  the  mesh.  For  example,  the 
user  could  partition  the  geometry  in  order  to  mesh  the  RVE  with  brick  elements. 
Otherwise,  the  add-on  will  automatically  mesh  the  RVE  with  tetrahedron  elements. 
The  add-on  then  inserts  double  nodes  along  the  delaminations  to  form  free  surfaces 
and  saves  the  mesh  and  surface  information. 


Figure  2.  Micrograph  of  CMC  showing  fracture  surfaces  [6]. 
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The  input  and  output  files  use  the  XML  (Extensible  Markup  Language)  format. 
The  format  is  supported  by  the  programming  languages  of  the  DDSHM  project 
(FORTRAN,  Mathematica™,  and  Python).  The  format  was  chosen  for  its  flexibility 
and  shallow  learning  curve.  Prior  knowledge  of  XML  is  not  required  to  use  the 
software.  A  Mathematica™  package  has  also  been  developed  to  create  the  input  files. 

The  resulting  hexagonal  array  RVE  included  six  damage  spaces  to  be 
implemented  in  the  UMAT:  one  defining  breaks  in  the  fibers,  three  defining  cracks 
in  the  matrix,  another  defining  separation  between  the  fiber  and  the  BN  coating  and 
lastly  a  damage  space  defining  the  separation  between  the  BN  coating  and  the 
matrix.  Six  unique  cracks  were  defined  in  the  RVE.  Each  crack,  or  damage  space, 
represents  a  different  mode  of  failure.  Figure  3  details  each  of  the  six  damage 
spaces  placed  in  the  RVE  and  the  associated  damage  modes. 


>.l  -  F  iber  Break  Cracks 
}2  -  Normal  Matrix  Cracks 
JJ  -  Matrix  Cracks 

>.4  -  Cracks  at  the  Interface  of  the  BN  and  Matrix 
L5  -  Cracks  at  the  Interface  of  the  F iber  and  BN 
>.6  -  Sleeve  Matrix  Cracks 


Figure  3.  Placed  cracks 
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The  resulting  UMAT  was  used  along  with  simple  FE  model  with  layered  elements 
oriented  to  represent  a  plies  in  the  0°  and  90°  directions  in  order  to  approximate  a 
woven  layer.  While  a  more  complex  woven  RVE  could  be  generated  [8],  for  this 
initial  study  the  simpler  hexagonal  array  was  utilized.  To  illustrate  the  functionality  of 
the  UMAT,  several  examples  are  provided. 

First,  with  reference  to  Eq.  (7),  note  that  q  is  dependent  upon  time,  where 

7(0  =  £-7o  (10) 

At0 

where  qo  is  the  assigned  viscosity  parameter,  Ato  is  an  arbitrary  time  step  and  At  is 
current  time  step  of  the  model.  This  allows  the  behavior  of  the  damage  to  be 
dependent  upon  on  the  rate  of  the  applied  load  and  allows  for  the  modeling  of 
creep.  Runs  with  a  linearly  increasing  load  for  three  different  values  of  q,  shown  in 
Figure  4,  illustrate  the  effect  of  this  parameter.  For  smaller  values  of  q,  the  crack 
requires  more  energy  to  open  fully.  Larger  values  of  q  require  less  load  to  open  fully 
and  open  to  the  full  extend  more  quickly.  Note  that  for  cases  shown  in  the  figure, 
smaller  values  of  q  result  in  slower  damage  growth  and  the  final  damage  state  is  not 
reached  over  the  range  of  load  analyzed. 

Creep  tests  were  also  modeled  to  demonstrate  the  ability  to  capture  crack  growth 
occurring  over  time  with  the  loading  held  constant.  Because  crack  growth  due  to 
temperature  is  an  important  consideration  for  CMCs,  a  model  with  ends  fully- 
constrained  was  used  to  demonstrate  crack  grow  under  applied  and  constant 

temperature  as  well  (note  that  for  simplicity  there  is  no  functional  dependence  of  G°r 
on  temperature  T).  In  both  cases,  the  load  was  ramped  for  ten  load  steps  and  then  held 
for  an  additional  ten  steps.  Normalized  crack  growth  for  the  fiber  crack  is  illustrated 
in  Figure  5  and  Figure  6  for  the  mechanical  and  thermal  loads,  respectively.  In  both 
cases,  the  crack  grows  as  expected  with  increasing  load  or  temperature  and  continues 
to  grow  with  time  as  the  load  or  temperature  is  held. 


Crack  Growth  for  the  FiberCrack  in  the  Element  Orientated 
inthex-Direction 


Applied  Load,  kips 


Figure  4.  Normalized  crack  growth  for  the  three  values  of  q 
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Figure  5.  Normalized  crack  growth  for  the  ramp  plus  hold  load  profile 


Crack  Growth  for  the  Fiber  Crack  in  the  Element  Orientated  in 
the  x- Direction 


0  100  200  300  400  500  600  700  800  900  1000  1100 


Applied  Temperature,  degF 


Figure  6.  Normalized  crack  growth  for  the  ramp  plus  hold  thermal  profile 


An  additional  analysis,  shown  in  Figure  7,  demonstrates  the  ability  to  apply  loads 
at  different  rates  and  shows  the  results  of  applying  load  at  two  different  rates.  The 
time  step,  dt,  of  1 .0  results  in  an  increased  rate  of  the  applied  load.  Thus,  the  crack 
opens  much  sooner  and  in  fewer  steps.  The  time  step  of  0.1  results  in  a  much  slower 
load  rate  and  crack  slowly  opens  over  a  number  of  time  steps.  The  models  eventually 
end  at  the  same  damage  state  but  the  time  it  takes  each  of  them  is  different  as  is  the 
applied  load  at  which  the  crack  is  fully  open. 
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Comparion  of  Crack  Growth  when  Applying  the  Load  at  Different 

Rates 


Applied  Load,  kips 


Figure  7.  Normalized  crack  growth  at  different  load  rates 


RESULTS  AND  DISCUSSION 


This  section  presents  a  sample  analysis  conducted  to  illustrate  the  implementation 
of  DDSHM  for  predicting  nonlinear  response  of  ceramic  matrix  composites.  The 
functionalities  demonstrated  above  allow  the  modeling  of  more  complex  behaviors 
of  materials.  However,  the  focus  of  the  initial  study  has  been  on  demonstrating 
integration  of  the  micromechanics  model  with  the  commercially  available  finite 
element  code  ABAQUS,  via  user-material  subroutines.  Toward  that  end  the  UMAT 
was  applied  to  a  simple  structural  element,  a  beam  in  4-point  bending  under 
monotonically  increasing  load. 

Before  moving  on  to  presentation  of  those  results,  it  is  important  to  note  that  the 
strain  increment  to  be  attempted  in  each  load  step  is  selected  by  the  main  ABAQUS 
program  based  upon  information  supplied  by  UMAT  for  the  previous  load  increments 
and  thus,  is  not  controlled  by  the  user.  The  stress  increment  returned  by  UMAT  is 
used  to  check  for  equilibrium  of  the  structure  at  the  global  level  (within  some 
numerical  tolerance).  It  is  easy  to  see  then  that  the  total  computational  time  in  the 
solution  phase  shown  in  Figure  1  depends  not  only  on  the  number  of  computations 
(related  to  N  for  DDSHM  and  mDOF  for  a  global-local  approach)  performed  at  the 
local  level  but  also  depends  on  the  stability  of  the  local  constitutive  routine.  Stability 
here  refers  to  the  ability  of  the  local  model  to  return  the  expected  updated  stress  vector 
and  stiffness  tensor  for  the  material  point  to  the  global  model  such  that  global 
equilibrium  is  obtained  in  the  fewest  iterations  possible.  An  example  of  an  unstable 
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local  model  is  one  where  the  stiffness  and  stress  residual  varies  in  a  discontinuous 
manner  for  small  values  of  state  variable  evolution.  Because  the  effective  moduli  in 
the  DDSHM  are  based  upon  a  piece-wise  interpolation  of  values  at  discrete  points  of 
crack  evolution,  smoothness  is  achieved  through  a  finite  difference  technique. 

As  outlined  above,  DDSHM  approach  allows  the  integration  of  a  nonlinear 
material  model  into  ABAQUS,  a  commercial  finite  element  code  for  the  structural 
analysis.  This  accomplished  by  inserting  the  nonlinear  material  as  a  user  defined 
material.  Using  the  RVEs  discussed  above,  ABAQUS  FE  models  were  developed  for 
the  4-Point  Bending  tests  shown  in  Appendix  R  of  [6],  The  fracture  data  available  for 
SiC  composites  indicates  a  fairly  wide  range  of  SERR  values  [9,  10].  Therefore  the 
Short  Beam  Shear  data  available  in  the  same  report  was  used  to  determine  the  values 
that  would  result  in  matrix  cracking  and/or  matrix/interface  cracking  at  the  point 
where  nonlinearity  becomes  apparent  in  the  test  data.  Another  advantage  of  these 
UMATs  is  that  layers  of  elements  can  be  defined  with  multiple  orientations  while 
referencing  the  same  UMAT  and  the  material  oriented  at  ±45  to  represent  the  woven 
material  in  the  test  specimens. 

The  UMAT  used  for  these  analyses  operates  by  storing  the  stiffness  of  an  RVE  for 
all  possible  combinations  of  crack  permutations.  As  the  elements  in  the  finite  element 
model  are  loaded  the  UMAT  calculates  the  strain  energy  release  rates  (SERRs)  for 
each  of  the  cracks.  When  the  SERR  reaches  a  defined  critical  SERR  for  a  particular 
the  crack,  that  crack  will  begin  to  open  and  the  stiffness  of  the  element  will  change  in 
the  finite  element  model  to  the  previously  determined  stiffness  which  matches  the 
current  state  of  damage  of  the  element.  Thus  for  the  next  step  of  the  loading,  the 
stiffness  of  the  damaged  elements  will  have  changed,  hence,  a  progressive  damage 
model.  The  crack  opening  for  each  crack  type  in  each  element  is  stored  as  a  state 
dependent  variable  and  can  be  plotted  in  post-processing.  This  allows  the  user  to  view 
damage  progression  by  crack  mode. 

The  load  vs.  displacement  curve  from  applying  the  derived  UMAT  to  the  4-point 
bending  model  is  compared  with  experimental  data  from  [6]  and  the  curve  for  the 
material  with  no  damage  in  Figure  8.  The  comparison  with  test  data  is  quite  good 
until  very  close  to  failure  (P  ~  80  lb).  Driving  the  model  past  convergence  issues  far 
into  progressive  failure  requires  loosening  convergence  parameters  too  much  to 
capture  the  material  response  however  this  effect  generally  provides  a  reasonable  and 
conservative  estimate  of  ultimate  load. 
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Deflection  (in.) 


Figure  8.  Comparison  of  4-pt  bend  test  data  for  S200-1  [6]  with  linear  material  and 

DDSHM  models 


Plotting  the  damage  parameters  allows  features  of  ceramic  matrix  composite 
failure  to  be  examined.  For  instance,  the  extent  of  transverse  matrix  cracks,  matrix 
sleeve  cracks,  matrix/interface  cracks  and  fiber  cracks  are  plotted  at  a  point  near 
specimen  failure  (P  ~  80  lb)  in  Figure  9.  Note  that  as  would  be  expected,  transverse 
matrix  cracking  is  more  extensive  than  anything  else  and  the  matrix/interphase  cracks 
have  grown  over  a  large  portion  of  the  specimen  while  fiber  cracking  is  far  more 
limited. 
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Figure  9.  Normalized  crack  extents  at  P  ~  80  lb  for  the  4-point  bending  FEM  using 

the  hexagonal  array  RYE  UMAT. 


CONCLUSIONS 

This  paper  has  presented  current  work  on  the  development  and  implementation 
of  representative  volume  elements  to  derive  the  three-dimensional  thermo¬ 
mechanical  constitutive  relations  for  ceramic  matrix  composites  with  multiple 
growing  cracks  using  the  discretized  damage  space  homogenization  method 
(DDSHM)  approach,  which  the  authors  believe  represents  the  most  promise  for 
advancing  the  state  of  the  art  in  durability  modeling.  The  constitutive  and  evolution 
equations  developed  from  the  RVE  have  been  integrated  and  implemented  in  a  general 
purpose  structural  analysis  code  (ABAQUS),  via  user-material  subroutines  (UMATs). 
Simple  FE  models  have  been  used  to  demonstrate  the  code’s  ability  to  capture  creep, 
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load  rate  and  thermal  effects  on  damage  growth  and  the  method  has  been  further 
demonstrated  using  an  ABAQUS  FE  model  of  a  4-point  bending  test  which  has 
provided  good  correlation  with  data  from  the  literature. 

While  this  establishes  in  principle  the  usefulness  of  this  approach  for  the  analysis 
of  ceramic  matrix  composites,  further  work  is  necessary.  In  particular,  continuing 
development  will  focus  on  coding  the  analysis  to  tie  a  thermo-oxidation  model  to  the 
crack  evolution  law  implemented  in  the  RYE. 
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